#读取序列化后文件并进行一定的判别处理
#x：产品名称
week_de<-function(x="product_1"){
  data<-list()
  #读取目标文件
  data$r$all <- read.csv(paste("DATA/",x,"_r_1",sep=""))
  data$p$all <- read.csv(paste("DATA/",x,"_p_1",sep=""))
  rownames(data$r$all) <- data$r$all[,1]
  rownames(data$p$all) <- data$p$all[,1]
  data$r$all <- as.matrix(data$r$all[,-1])
  data$p$all <- as.matrix(data$p$all[,-1])
  #空数据统计标记
  data$pr$zero<-(data$p$all!=0|data$r$all!=0)%*%matrix(1,nrow = (ncol(data$p$all)))
  data$pr$zero<- cut(data$pr$zero,breaks = c(-Inf,ncol(data$p$all)/3,ncol(data$p$all)*2/3,Inf),labels = c("l","m","h"))
  table(rownames(data$p$all)==rownames(data$r$all))
  #店铺销售均值
  data$pr$median<-apply(data$p$all-data$r$all,1,median)
  data$pr$median<-cut(data$pr$median,breaks = c(-Inf,5,20,Inf),labels = c("l","m","h"))
  #全销售数据统计标记
  data$pr$full<-apply((data$p$all-data$r$all)<=0,1,sum)
  data$pr$full<-data$pr$full/ncol(data$p$all)
  data$pr$full<-cut(data$pr$full,breaks = c(-Inf,0.2,Inf),labels = c("h","l"))
  data
}
#读取环境中的序列数据并进行治理与数据选择(日为单位，周为基单位自动治理)
#x：序列集
#n：序列ID
#start:序列起始周期
#end:序列终止日期
#k:序列周期基数
#返回序列拆分结果以及序列剔除星期数
n=6
x=data
lgts<-function(x,n,start=1,end=37,k=2){
  #计算数据序列总长度
  date_p<-strptime(substr(colnames(x),start = 2,stop=11),format="%Y.%m.%d")
  date_f<-strftime(date_p,format = "%Y-%m-%d")
  #最大日期与最小日期差天数
  ltime<-difftime(min(date_p),max(date_p),units="days")
  ts<-as.data.frame(cbind(date_f,x[n,]))
  colnames(ts)<-c("时间","频数")
  #制作日期序列
  dates<-as.data.frame(strftime(as.Date.numeric(x = 1:(1-ltime),origin = strptime("2012-07-01",format="%Y-%m-%d")),format = "%Y-%m-%d"))
  colnames(dates)<-"时间"
  ts<-join(x=dates,y=ts,by = "时间",type = "full",match = "all")
  #这里是日期截取start=x时截取(x,1)，end=x时，截取(x,7)。
  ts<-ts[((start-1)*7+1):(end*7),]
  flag<-table(weekdays.Date(strptime(ts[is.na(ts$频数),]$时间,format="%Y-%m-%d")))
  rm<-c(names(flag)[flag>(0.1*length(ts$频数))])
  if(length(rm)!=0){
    ts<-ts[!(weekdays.Date(strptime(ts$时间,format="%Y-%m-%d"))%in%rm),]
  }
  ts$频数<-as.numeric(as.character(ts$频数))
  ts$频数[is.na(ts$频数)]<-0
  #构建时间序列
  cyc<-(7-length(rm))
  timed<-ts(ts$频数,start=c(floor((start-1)/k)+1,((start-1)%%k)*cyc+1),end=c(floor((end-1)/k)+1,((end-1)%%k+1)*cyc),frequency =k*cyc)
  timed_t<-decompose(timed)
  plot(timed_t)
  list(timed_t,rm)
}
#绘制定阶图的函数(需要以周为单位)
#ts：想要绘制的时间序列
#lag.max:时间序列观测长度
#lgval：自动定阶系数，未用到
#dif：周期长度
#sen：是否显示周期定阶图
#diff：差分阶数（目前只实现一阶）

plts<-function(ts=ts,lag.max=60,lgval=0.15,dif=8-length(ts[[2]]),sen=TRUE,diff=0){
  if(diff!=0){
    tmp<-ts[[1]]$x[(1+diff):length(ts[[1]]$x)]-ts[[1]]$x[1:(length(ts[[1]]$x)-diff)]
  }else{
    tmp<-ts[[1]]$x
  }
  #自相关系数检验
  pic1<-acf(tmp, lag.max, type=c("correlation"),plot=FALSE)
  #篇相关系数检验
  pic2<-acf(tmp, lag.max, type=c("partial"),plot=FALSE)
  #进行季节性差分
  dtime<-tmp[(dif+1):length(tmp)]-tmp[1:((length(tmp)-dif))]
  #ts.plot(dtime)
  #自相关系数检验
  pic3<-acf(dtime, lag.max, type=c("correlation"),plot=FALSE)
  #偏相关系数检验
  pic4<-acf(dtime, lag.max, type=c("partial"),plot=FALSE)
  if(sen==TRUE){
    pic<-as.data.frame(rbind(cbind(pic3[[1]][-1],"SACF",1:lag.max),cbind(pic4[[1]],"SPACF",1:lag.max)),stringsAsFactors = FALSE)
  }else{
    pic<-as.data.frame(rbind(cbind(pic1[[1]][-1],"ACF",1:lag.max),cbind(pic2[[1]],"PACF",1:lag.max)),stringsAsFactors = FALSE)
  }
  colnames(pic)<-c("val","class","lag")
  pic$val<-as.numeric(pic$val)
  pic$lag<-as.numeric(pic$lag)
  #绘制相关图
  p<-ggplot(as.data.frame(pic),aes(x =lag,y=val))+
    geom_bar(stat="identity",fill="#00AEAE",colour="black")+
    facet_wrap(~class,nrow=2)+ylim(c(-1,1))+
    geom_hline(yintercept=c(-lgval,lgval))+
    labs(title = "ARIMA模型主要定阶指标")+
    theme(plot.title = element_text(size = rel(2)))+
    labs(x = "滞后数", y = "相关系数")
  res<-c(table(abs(pic1$acf[-1])>lgval)[2],
         table(abs(pic2$acf)>lgval)[2],
         table(abs(pic3$acf[-1])>lgval)[2],
         table(abs(pic4$acf)>lgval)[2])
  res[res>lgval*lag.max]<-0
  list(res,p)
}

#预测模型
#data：序列数据
#id:序列编号
#start:序列起始位置
#end：序列终止周期
#k：序列周期阶数
id=4

##应重写短周期转长周期函数
myaim<-function(data,id=6,start=17,end=28,k=4){
  #通过lgts，获取训练序列
  ts <- lgts(data,n=id,start,end,k)
  #获取序列周期
  cyc <- (7-length(ts[[2]]))
  #计算序列趋势
  y<-lgemc(data=ts[[1]]$x,cyc)
  #计算序列长趋势，与去除长趋势后的残差
  yfit_s4<-loemc(y,cyc)
  plot(as.numeric(ts[[1]]$x),type="l")
  lines(yfit_s4,col="red")
  tsx<-ts[[1]]$x-yfit_s4
  tsx<-ts(tsx,start=c(floor((start-1)/k)+1,((start-1)%%k)*cyc+1),end=c(floor((end-1)/k)+1,((end-1)%%k+1)*cyc),frequency = cyc*k)
  tsxp<-decompose(tsx)
  plot(tsxp)
  rd <<- k*cyc
  ARIMA<-arima(tsx,order = c(rd/k,0,0),seasonal=list(order=c(2,1,1),period=rd),method = "CSS")
  #ARIMA<-arima(ts[[1]]$x,order = c(2,0,0),seasonal=list(order=c(2,1,1),period=rd),method = "CSS")
  pre=predict(ARIMA,n.ahead=rd)
  pre$pred<-pre$pred+yfit_s4[1:rd]
  U=pre$pred + 1.96*pre$se 
  L=pre$pred - 1.96*pre$se
  tsp <- lgts(data,n=id,start,end+4,k)
  tsp_x<-tsp[[1]]$x
  #length(tsp_x)
  tsp_x<-as.data.frame(cbind((((start-1)*cyc+1):(((end+4)*cyc)))/rd,tsp_x))
  colnames(tsp_x)<-c("时间","销售量")
  LU<-as.data.frame(cbind((((end)*cyc+1):(((end+4)*cyc)))/rd,L,U,pre$pred))
  colnames(LU)<-c("时间","L","U","pred")
  pp<-ggplot(data=tsp_x,aes(x=时间,y=销售量))+
    geom_line()+geom_line(data=LU,aes(x=时间,y=pred),colour="red")+
    geom_polygon(aes(x=c(时间,时间[rd:1]),y=c(L,U[rd:1])),data=LU,alpha=0.3)+
    xlim(c(end/k-5,end/k+1))+
    labs(title = paste("店铺",rownames(data)[id]))+
    theme(plot.title = element_text(size = rel(2)))+
    geom_vline(xintercept =((end+(1/cyc))/k), colour="green", linetype = "longdash",size=1)
  list(pp,ts[[1]]$x,ARIMA)
}
#拟合长趋势值
loemc <- function(y,cyc){
  if(which.max(y[,2])!=length(y[,2])){
    ymax <- y[c((which.max(y[,2])+1):length(y[,2]),1:which.max(y[,2])),]
  }else{
    ymax <- y
  }
  #构建拟合模型
  model <- lm(ymax[,2] ~ poly (1:length(ymax[,2]),7))
  yfit<-cbind(ymax,fitted(model))
  yfit_sort<-yfit[order(yfit[,1]),]
  yfit_s4<-matrix(matrix(1,nrow=cyc)%*%matrix(yfit_sort[,3],nrow=1),ncol=1)
  yfit_s4
}
#计算移动平均趋势值
lgemc <- function(data,cyc=4){
  data<-as.data.frame(cbind(floor(((1:length(data))-1)/cyc),data))
  colnames(data)<-c("周数","数量")
  res<-aggregate(data$数量,by=list(data$周数),mean)
  colnames(res)<-c("周数","数量")
  res
}
plot(ts)
#周刊预测
myaim_wk<-function(data=tmp,id=2,start=12,end=23,k=4){
  names<-rownames(data)[id]
  data<-data[id,data[id,]!=0]
  data_t<-data[((start-1)*k+1):(end*k)]
  wts <- ts(data_t,start = start,end=c(end,k),frequency = k)
  wts[wts>(3*sd(wts)+mean(wts))]<-mean(wts)
  ts <- decompose(wts)
  #plot(ts)
  
  #model <- lm(ymax[,2] ~ poly (1:length(ymax[,2]),7))
  #yfit<-cbind(ymax,fitted(model))
  #yfit_sort<-yfit[order(yfit[,1]),]
  #yfit_s4<-matrix(matrix(1,nrow=cyc)%*%matrix(yfit_sort[,3],nrow=1),ncol=1)
  #tsx<-ts[[1]]$x-yfit_s4
  #tsx<-ts(tsx,start = 16,end=c(28,cyc*k),frequency = cyc*k)
  rd <<- k
  ARIMA<-arima(ts$x,order = c(0,0,0),seasonal=list(order=c(1,1,1),period=rd),method = "CSS")
  #ARIMA<-arima(ts[[1]]$x,order = c(2,0,0),seasonal=list(order=c(2,1,1),period=rd),method = "CSS")
  pre=predict(ARIMA,n.ahead=rd)
  # pre$pred<-pre$pred+yfit_s4[1:rd]
  U=pre$pred + 1.96*pre$se 
  L=pre$pred - 1.96*pre$se
  tsp <- ts(data[((start-1)*k+1):((end+1)*k)],start = start,end=c(end+1,k),frequency = k)
  tsp_x<-tsp
  tsp_x[tsp_x>(3*sd(tsp_x)+mean(tsp_x))]<-mean(tsp_x)
  tsp_x<-as.data.frame(cbind(((start*rd):(((end+2)*rd-1)))/rd,tsp_x))
  colnames(tsp_x)<-c("时间","销售量")
  LU<-as.data.frame(cbind((((end+1)*rd):(((end+2)*rd-1)))/rd,L,U,pre$pred))
  colnames(LU)<-c("时间","L","U","pred")
  pl<-ggplot(data=tsp_x,aes(x=时间,y=销售量))+
    geom_line()+geom_line(data=LU,aes(x=时间,y=pred),colour="red")+
    geom_polygon(aes(x=c(时间,时间[(rd):1]),y=c(L,U[(rd):1])),data=LU,alpha=0.3)+
    xlim(c(start,end+2))+
    labs(title = paste("店铺",names))+
    theme(plot.title = element_text(size = rel(2)))+
    geom_vline(xintercept =(end+1), colour="green", linetype = "longdash",size=1)+
    geom_hline(xintercept = 0,colour="red")
  list(ARIMA,pl,ts$x)
}



